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Abstract 

We study the sojourn time in a queueing system with a single exponential server, serving a Poisson stream 
of customers in order of arrival. Service is provided at low or high rate, which can be adapted at exponential 
inspection times. When the number of customers in the system is above a given threshold, the service rate is 
upgraded to the high rate, and otherwise, it is downgraded to the low rate. The state dependent changes in 
the service rate make the analysis of the sojourn time a challenging problem, since the sojourn time now also 
depends on future arrivals. We determine the Laplace transform of the stationary sojourn time and describe 
a procedure to compute all moments as well. First we analyze the special case of continuous inspection, 
where the service rate immediately changes once the threshold is crossed. Then we extend the analysis to 
random inspection times. This extension requires the development of a new methodological tool, that is 
matrix generating functions. The power of this tool is that it can also be used to analyze generalizations to 
phase-type services and inspection times. 

Key words. Sojourn time distribution. Matrix generating function, Adaptable service speed 
AMS subject classifications. 60K25,60K37,60D05. 


1 Introduction 

We consider a single-server queueing system, where customers arrive according to a Poisson stream with rate 
A and receive service in order of arrival. The service requirements are exponential with mean 1. The rate of 
the server can be either /tq, or /ri and this service rate can be adapted at random inspection times that occur 
according to a Poisson stream with rate 7 . For convenience, we think of as the fastest rate, i.e. /ri > even 
if under the stability conditions this assumption may be removed. When the number of customers in the system 
is above the threshold K, the service rate is upgraded to the high rate /ri, and otherwise, it is downgraded to 
the low rate /tq- An important performance measure is the sojourn time. In this paper we aim to determine 
its stationary distribution. This is a challenging problem, since due to adaptable service rate, the sojourn time 
does not only depend on the state seen at arrival, but it also depends on future arrivals. 

There is a considerable literature on the analysis of single server queueing systems with variable service rates, 
see e.g. [DlEllIlIIIllIl]. Those studies often assume that the service rates can be continuously adapted based 
on the queue content, and focus on the calculation of the steady-state workload distribution. An exponential 
multi-server system is considered in m, with the feature that a reserved block of servers can be switched on 
(which takes an exponential switch-on time) when the number of customers in the system exceeds a certain 
threshold, and this block is immediately switched off when the number drops below another threshold. The 
emphasis in |15j is on the trade-off between the mean sojourn time and operating costs of the servers. In [3], the 
stationary distribution of the workload is determined for an MjGjl queue, where the service rate can not be 
continuously adapted, but only right after customer arrivals. In the literature, systems with adaptable service 
speed at inspection times have already been analyzed, we refer the reader to [1[5] and the references therein. 

The model with inspection rate 7 < 00 can be handled by considereing a two-stage birth-death process. This 
kind of model usually shows up in the analysis of retrial queues, where the state of the system has to keep track 
of the size of the retrial orbit. We refer the reader to the survey [1^. In [9], the number of retrials of a generic 
customer is analyzed, which is a quantity directly related to the sojourn time and which depends on future 
arrivals to the system. Paper starts the analysis with a matrix equation that is similar to the one appearing 
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in Section O but it is able to reduce this equation to a scalar one by exploiting the fact that the retrial system 
has no buffer and only one server. Multi-channel systems are much more complicated to analyze and very little 
results are available about the sojourn time. Generally, what makes retrial systems more complicated than the 
one analyzed here, is the property that the rate at which retrial customers arrive to the system is proportional 
to the size of the orbit. This phenomenon is not appearing in our system, which is one of the reasons why our 
analysis is feasible. 

As mentioned above, the focus in the current paper is on the sojourn time, not on the workload or number 
of customers in the system. In Section [5J we first consider the special case of continuous inspection (so 7 = 00 ), 
where the service rate immediately changes once the threshold K is crossed. This assumption simplifies the 
model, though it still contains the complication that the sojourn time depends on future arrivals. For the case 
of continuous inspection, we determine the Laplace transform of the stationary sojourn time and describe a 
procedure to compute all moments as well. The computation of the Laplace transform requires a recursive 
scheme and for the case 7 < c» the Laplace transform can be expressed in terms of matrix functions that can 
be computed as solutions of a linear matrix system. 

Then, in Section [31 we proceed by extending the analysis to random inspection times occurring according a 
Poisson stream with rate 7 < 00 . This extension, however, is not straightforward, but it requires the development 
of a new methodological tool, that is matrix generating functions. By employing this tool we are able to find 
an expression for the Laplace transform of the stationary sojourn time, involving finitely many terms which can 
be recursively calculated. The analytical results are illustrated by numerical examples. 

2 Model with continuous inspection 

In this section we first consider the special case of continuous inspection, so 7 = 00 . This implies that whenever 
the number of customers in the system exceeds the threshold K > 0, the rate of the server is immediately 
upgraded from the low rate /xg to the high rate fii > fj-Q. As soon as the number of customers in the system 
becomes less or equal to K, the rate of the server is reduced to the low rate /xg again. 

Denoting by Q(t), the number of customers in the system at time t > 0, we have that the process is a 
continuous time Markov chain, the transition diagram of which is depicted in Figure [I] 



Figure 1: Transition diagram for continuous ispection model 


Denoting by Q* the stationary number of customers in the system, we have that its distribution is given by 


TTn = P(Q* =n) = 


(A//xg)” 7 rg for n < AT 


(/Xi//Xg)^ TTg for n > a: 

and under the stability assumption //i > A, the value of ttq is given by 

/ K 


TTo = y](A//io)" 


A 


\n—0 


/Xi - A 


(Vmo) 


K 


( 1 ) 


( 2 ) 


We aim to compute the distribution of the sojourn time of a typical customer that arrives to the system in 
stationary regime. Note that, in order to do this, we can not use Little’s distributional law |13| . since future 
arrivals may affect the sojourn times of the customers already present in the system by inducing a change in 
the service rate. 

As shown in Figure |21 we identify a tagged customer in the queue by a pair of numbers where n 

stands for the position of the tagged customer in the queue, and where m denotes the number of customers 
behind him. We denote the sojourn time of this customer {n,m) by S(n,m). The stationary sojourn time is 
denoted by S*. 



customer 



Figure 2: Tagged customer {n,m) at position n in the queue, with m customers behind him. 
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For the Laplace transforms ip{s) = E[e ] and 'il}{n,m,s) = E[e the following relation holds by 

virtue of PASTA [T7] , 

OO 

i’i.s) = '^^{n + l,Q,s)TTn. (3) 

n=0 


Hence, to compute the Laplace transform of the stationary sojourn time S*, we need to compute the transforms 
'(/'(n, 0 , s) for each n > 0 . 

By using next-event analysis we have, for n > 0, 


S{n, m) = < 




X+ 


A-t-/2] 


S{n — 1, in) 

w.p. 

Ao/ (ao + A) 

S{n, TO -1- I) 

w.p. 

A/(ao + A) 

S{n — 1, to) 

w.p. 

Ai /(ai + A) 

S{n, m + 1) 

w.p. 

A/ {fJ-i + A) 


, as n -|- TO < AT 


, as n -I- TO > AT 


(4) 


where X denotes an independent exponential random variable with rate 1, and S'(0 ,to) = 0. By Laplace 
transforming the relations ®, we get, for n > 0 , 


ipin, TO, s) 


f^i{n+m>K} fp{n -l,m,s) + X ip{n, m + l,s) 
^ T f^l{n+m>K} T ^ 


(5) 


with boundary conditions, '!/'( 0 ,to, s) = 1 , for all to > 0 , and where we used the indicator function 1 {A} = 1 if 
A is true and 0 otherwise. 

When TO > A', it follows that for any n > 0, the server will work at high speed during the whole sojourn 
time of the (n,TO)-tagged customer. Hence S{n,m) is Erlang distributed with parameters n and fii, and thus 
its Laplace transform is equal to 


ip{n,m,s) = + s))'^ as n > 0 and TO > AT. ( 6 ) 

The above equation is also valid for n = 0. 

Using expression in ®, the Laplace transforms ip(n, to , s) for m < K can be recursively computed in n, 
as the following lemma shows. The proof of the lemma is deferred to Appendix lAl 

Lemma 2.1. By defining 


as{k) =^o/(A + f^o + s)l{k < K} + fii/ {X + fii + s)l{k > K} ] 
&s(fc) =A/(A -|- fiQ s)l{/c < AT} -|- A/(A -I- -l-s)l{fc>Ar}, 

and Bs{k, 0 ) = 1 and Bs(k, h + 1) = Bfik^ h) bs{k + h) for k,h > 0, we have 


tpin, TO, s) =Bs(n + m,K — to) 
K-l 


Ai 


g-i + s ^ 

+ as{n -I- k) Bs{n + m,k — to) ip(ji — 1 , k, s ), 


(7) 


k—m 


for n > 0 and 0 < m < K. 

Remark. It can be easily shown that the value of Bfik, h) can be explicitly computed by the following formula 


Bs{k,h) = 


A 


s -l- A -l- fii 


A /s-l-A-l-/ii\ 
/ V '5 + A-l-/io/ 


hA{K-k+l) + 


( 8 ) 


with a Ab = min{o, b} and (a)"'" = max{a, 0}. 

Relation (jS]) and Lemma l2 .1 1 allow us to compute s) for any to, n > 0. However, to calculate ip{s) in 

m we still need to compute an infinite number of terms. To overcome this issue we take advantage of the fact 
that, above the threshold AT, the transition diagram is invariant towards the right, similarly to the standard 
M/M/1 queue. To use this invariant property we introduce the following marginal z-transform 


(j>{z, TO, s) = + h + l,m^s) , 

h=0 


(9) 


valid for \z\ < 1. In the following we show how to compute, in finitely many steps, the function (j){z,m, s). We 
use it to calculate the infinite sum in m and then obtain a formula to compute the Laplace transform of the 
sojourn time as given in Proposition 12.21 
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By writing ([S]) for n = K + h + 1, multiplying by and summing over all /i > 0, the following recursive 
equation holds 

Hi'tl){K,m,s) + l,s) 

</)(z,m,s) = 


A + /ii(l — z) + s A + /ri(l — z) + s 


( 10 ) 


As boundary value we have 

CX3 

4>{z,K,s) 


h=0 


Ml 


K+h+1 


Ml 


Ml + s, 

K+l 


= 


Ml 


Ml + s 


K-\-l oo 


E 

h=0 


Ml 

Ml + s 


Ml + s 


^Mi + sJ ni{l — z) + s ' 
from which the values of 4){z^ m, s) can be recursively computed for m = AT — 1 ,..., 0, yielding 


( 11 ) 


K-l-r 


A 


h^O 


A + Mi(l — z) + s 


K-r. 


(t>{z,K,s). 


( 12 ) 


In particular we can compute, in finitely many steps, the value of (^(z, 0, s). 

Knowing the value of (j){z,Q,s), expression (jS]) can be finally computed as summarized in the following 
proposition. 


Proposition 2.2. The Laplace transform of S* can be eomputed in the form 


K-l 


V'(s) =7^0 


h—0 _ 


A 


— ) V'(l7 + 1) Oj s) + ( — ) ( — 

Mo/ VMo/ Vmi 

K 


A 


K 


A 


Ml 


Ml + s 


h+l 


ip{K,h,s) 




Ml 


Ml 


Ml - 

Ml/ + fii — X + s 


(13) 


Proof. The result follows from ([3|) by splitting the sum in a finite, n < K, and infinite part. 


■^-1 / y \ " f 

l/’(s) =710 E — •!/'(n + l,0,s) + 7ro — p,i,Q,s) 

^0 XToJ 


(14) 


For the last term we use m and m to get 

K — l — m 


i\ —1 —m / ^ \ 

(/)(A/Mi,m,s) = V — ) 


+ 


Ml 


h / \ h-\-l 

Ml 

Ml + s 


A 


Ml — A + s \Mi 


m + h,s) 
Ml 


Ml + s 


and the result follows by rearranging terms. 


(15) 

□ 


The terms appearing in equation o have the following nice probabilistic interpretation. 

• With probability nh = 7ro(A/Mo)^, h < K, the tagged user enters a system with h customers and experi¬ 
ences a sojourn time, the Laplace transform of which is 'ijj{h + 1, 0, s). 

• With probability ttk+h = 7ro(A/Mo)^(A/Mi)^) 0 < M < AT, he finds K + h customers waiting. We slightly 
modify the system and assume that the tagged customer overtakes M -|- 1 customers and occupies position 
AT instead of the last one in the queue. In addition, the server first serves the last M -|- I customers. Since 
the speed of the server depends on the number of customers waiting and not on their specific order of 
service, the first M -|- I services will be at speed mi taking an Erlang time with parameters M -|- I and mi 
to complete. What is left is the service time of the tagged customer, the Laplace transform of which is 
if{K,h,s). 

• With probability n> 2 K = 7ro(A/Mo)^(A/Mi)^(Mi/(Mi~'^)) the tagged customer finds at least 2K customers 
waiting. As before, he is going to occupy position AT, the sojourn time of which is Erlang distributed with 
parameters AT and mi- The number of customers he has overtaken is at least AT, and the time it takes to 
complete their services is the sum of AT exponential random variables with parameter mi plus a generic 
sojourn time of an M/M/I queue having mi as service speed. This last quantity is exponentially distributed 
with parameter mi — A. 
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Remark. From the Laplace transform of the stationary sojourn time given in ([3]) an explicit expression for the 
distribution can be obtained. Indeed, the inverse transformation is straightforward as the Laplace transform is 
a rational polynomial, the poles of which are all located on the real axis. To be more precise, the locations of 
the poles belong to the set 

A = {—(A + ^i), —(A + /xo), —^ 1 , —(/xi — A)} , 

implying that the density function is given by a linear combination of terms for a & A and k = 

0 ,1 ,..., mult(a) — 1, where mult(a) denotes the multiplicity of pole a. 

Remark. If we let K ^ oo in o , we recover ©• For any n > 0, ipln + l, 0, s) becomes the Laplace transform of 
an Erlang distribution with parameters n+1 and /xq, and ip{s) reduces to the Laplace transform of an exponential 
distribution with parameter /xq — A, that is the distribution of the sojourn time of a classical M/M/1 queue 
with service rate /xq. 

Remark. If IF = 0, only the last term in m is different from zero. Substituting tto = (/xi — A)//xi, given in 
we get that x/’(s) is the Laplace transform of an exponential distribution with parameter /xi — A, that is the 
distribution of the sojourn time of a classical M/M/1 queue with service rate ^i. 


2.1 First moment calculation 


As mentioned in Remark [31 it is possible to compute the distribution of the sojourn time, but it is easier to 
compute the moments by using the relation E[S'^] = (—(0+). In this section we show how to compute 
the first moment. However, by taking higher order derivatives of the Laplace transform, recursive expressions 
can be obtained to compute all moments. 

Let V = E[5'] and Vn^m = E[5'(n,m)]. With n > 0, from ([ 6 ]) we have for m> K, Vn,m = n/^i, and using 
d?]), for 0 < 77 X < AT, 


Vn m =—^ 0 + in + m,K — m) — BA (n + m, K — m) 

Ml 

K-l 

+ ^2 ao+ + fc) Bq+ (n + m, fc — to) 

k—m 

— aQ+ (rx + k) Bq+ {n + m,k ~ "m-) 

— ao+ (n + k) Bq+ {n + m,k ~'m-)^ ) 

where 

a'o+ik) =-/xo/(A + /io)^I{fc < K} - fj,i/{X + fj.ifl{k > K} ; 
b'^+ik) = - A/(A + fio)A{k <K}- A/(A + fii)A{k > K} , 

and B'^^i^k, 0) = 0 and B'^^{k, h-\-l) = B'^^{k, h) bo{k + h) + Bo{k, h) b'^^{k + h). 

The following algorithm shows how to recursively compute Vn,m for 0 < m < K: 

for i=l to n do 

for j=l to K-m do 

I compute Vi,K-j 

end 

end 

Algorithm 1: Computing Vn^m for n > 0 and 0 < m < K 
Finally, by applying Proposition [23 we get 



A^ 


I^h+1,0 H- fF 

Mo 



+ TTq 



/ 2K 

VMi/ V(Mi-A) 


h 


(y^K,h + 


h + l \ 

Ml J 


Ml 

(A-Ml) 



(16) 


(17) 


3 Model with inspection times 

In this section we analyze the system where inspection times occur according to a Poisson stream with rate 
7 < oo. So, in this case, there is no continuous inspection and adaptation of the service rate is delayed (with 
an exponential time) when the number of customers in the system crosses the threshold K. If at an inspection 
time the system is found congested with more than K customers, the service rate is immediately set to the fast 
rate mi and otherwise, if at most K customers are present, the service rate is set to the low rate mo- 
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Now we need to include the service rate in the state description of the system, resulting in the Markov chain 
shown in Figure [3] Note that for any number of customers in the system, the service rate can be high and low. 




Figure 3: Transition diagram for exponential inspection times 


Denoting by M. the stationary random service rate, let ttok = IP(Ad = /ro,Q* = n) and 7ri„ = = 

Q* = be the stationary probabilities to find n customers in the system with the server working at rate 
and fii, respectively. In what follows, the quantity 7r„ denotes the column vector with components (TTon, Trin)'*', 
where (•)^ is the transposition operator. The stationary distribution satisfies the balance equations 


— HiTTq + M TTl =0 

A7r„_i — H 2 TTn + M 7r„_|_i = 0 1 < n < K ; 

A 7r„_i - H 3 TTn + M TTn+l =0 n> K , 


where the transition matrices are defined by 


M = 


Mo 0 \ 

0 ) 


A = 


A 0 \ 

0 A j ’ 


and Hi = A + r 2 , H 2 = M + A + r 2 and = M + A + Fa, where 


r2 = 





(18) 


From the theory on quasi-birth-death processes [161 [14], we conclude that for n> K, the stationary proba¬ 
bility vector TTn can be written in the form 


TTK+h = R’" T^K, h>0 , 


(19) 


where the matrix R is the minimal non-negative solution of the matrix equation 

A-H3R + MR^ = 0. ( 20 ) 

Using the probabilistic interpretation of i?, or by solving the matrix equation (EOl), it follows that R is of 
triangular form and in particular, it is equal to 


R = 


Roo 0 

7 -Rqo a 

Ml 1 —-Roo Ml 


with Rqo = 


/ Mo+7+A 
I 2/j.o 


) 


2 


Mo * 


The value of ttk can be computed by the normalizing equation 


K-l 

y] e TTfe -I- e (J - R)~^ttk = 1 , 

/c=0 


( 21 ) 


( 22 ) 


with e the all-one (row) vector. 

By PASTA, as in ([3]), the Laplace transform of the stationary sojourn time is given by 

00 

= '^1p{n+ 1,0, s) TTn, (23) 

n—0 

where ip{ri, m, s) denotes the row vector {ipoin, m, s),ipi{n, m, s)), with ipiin, m, s) being the Laplace transform 
of the sojourn time Si{n,m) of a tagged customer who is at position (n,m) and the service rate is 1 = 0,1. 
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By using next-event analysis, we get the following recursive equations for the sojourn times, Si{n, m), i = 0,1, 
n > 0 


^ ( Si{n-l,m) ^.Y>. +Hi+l) 

Si{n,m) =- ---h < Si{n,m+1) w.p. X/{X + + j) (24) 

[ Si{n+m>K}{n,m) w.p. 7/(A-hMi+ 7 ) 

where X denotes an independent exponential random variable with rate 1, and ^^(O, m) = 0. Taking the Laplace 
transform of (|M1) yields the equation 

^{n,m,s) {H{s) - = '0(n - l,m,s) M + tj}{n,m + l,s) A (25) 


for n > 0, where 

H{s) = {^ + s)I + A + M , ro=l^l ly = 

and '0(Oj ■s) = C) with e the all-one (row) vector and / the identity matrix. 

Similarly to ([HI), when m > K and for any n > 0, we have that whenever an inspection occurs, the service 
rate is set and kept at the value /ii till the end of the sojourn time of the tagged customer. This implies that 
ip{n, m + 1, s) = 'ip{n, m, s) ior m > K and n > 0, and substitution in (1251) yields 

Tpin, m,s) = e T'^{s) as n > 0 and m> K , (26) 

with T(s) = M {H{s) - Ti - A)-h 

Equations (^Hj) and (unD allow us to get the values of ^(n, m,s) for any n,m > 0. However to compute 
expression (ESI) in finitely many steps, we still need to find a way to handle the infinite sum. So far, the analysis 
proceeds as in Section El and thus the next step would be to introduce the marginal 2 -transforms corresponding 
to (El), that is, (j)i{z,m,s) = + h + l,m,s)z^. However, this approach immediately fails, since the 

stationary probability distribution dm calls for a matrix generalization. The main contribution of this work is 
to provide this generalization by the introduction of the following matrix generating function, 

CO 

4>{Z, m, s) = E i}{K+ h+l,m,s)Z^ , (27) 

h—0 

where Z is any matrix with eigenvalues contained in the open unit disk of the complex plane. 

Remark. Since the absolute value of the Laplace transform 'ijji{n,m, s) is less or equal to one, the assumption 
on the eigenvalues of Z implies that the matrix generating function 4>{Z, m, s) is well defined. 

Let us rewrite expression (1251) for n > K in the alternative form, 

tj}{n, TO, s) = tfj{n — 1, TO, s) Tm{s) + fj{n, m+l,s) Ta(s) (28) 


with Ta{s) = A {H{s) — Ti) A G {A, M}. Multiplying expression (1^ on the right by for n = iL -I- /i-I- 1, 
and then summing over h > 0 and using that T Z^ T~^ = {T Z T~^)^ , we get a recursive equation for (/>(Z, to, s), 

(l){Z, TO, s) =i’{K, TO, s) Tm{s) + (j){TM{s) Z T^^{s), to, s) Tm{s) Z 

+ 4>{Ta{s) ZT]f\s),m+l,s)TA{s) . (29) 

The main difference between equations dm and (1291) is that in the latter we loose the commutative property 
of the product and the functions (j) need to be evaluated for different values of their arguments. The boundary 
condition is obtained from dm, 


f>{Z,K,s)=eT^+\s) ij2T^{s)Z^ = eT^+\s) S{Z,I,T{s)), 


\h^0 


with I being the identity matrix and where we employed the definition, 


S{Z,A,B) AZ’^ . 

h—0 


(30) 


(31) 


The matrix S{Z,A,B) is well defined for any matrix Z,A,B with Z and B having all eigenvalues inside the 
closed and open disk, respectively (so that the series converges). Note that T{s) in (^01) has all eigenvalues 
inside the open unit disk. The matrix S{Z,A,B) can be computed as the solution of a matrix equation as 
shown in the following lemma. The proof is deferred to the appendix. 
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Lemma 3.1. Let Z, A and B be three matriees with Z and B having all eigenvalues in the elosed and open 
disk, respectively, then the matrix function S = S{Z,A,B) is the unigue solution of the following matrix system, 

S-BSZ = A. (32) 

The next proposition shows that, in order to compute the Laplace transform of the stationary sojourn time 
S* in terms of a finite number of addends, only the value of 4>{R, 0, s) is needed. 

Proposition 3.2. The Laplace transform of S* can be eomputed in the form 


K-l 

if{s) = ^ ij)(n + 1, 0, s) TTn + 4>{R, 0, s) ttk ■ (33) 

71 — 0 

Proof. The result follows from (12311 by splitting the sum in a finite, n < K, and infinite part. For the latter 
part, we express ttk+h as in dm ) for h > 0 , and apply definition ( 1221 ). □ 

The computation of (fiR, 0, s) requires some additional machinery with respect to the one developed in 
Section [2] for the scalar case. Before giving the statement of the main result we need the following technical 
lemma, the proof of which is deferred to the appendix. The lemma states that the infinite sum of matrices 
appearing at the left-hand side of pip can be recognized as a matrix function S, which can be computed from 
the matrix system (IS21)- 

Lemma 3.3. Let Z, A and B be three matrices with Z and B having all eigenvalues in the elosed and open 
disk, respectively, and let Ti and T 2 be invertible matrices with Ti having all eigenvalues in the open disk, then 
the following relation holds, 

CO 

S{T 2 T’f Z tY B) T 2 Z'^ = S{Z, A T 2 S{Z, I,Ti),B). (34) 

h—0 

The following result allows us to compute the value of 4>{R, m, s) in finitely many steps. 

Theorem 3.4. The values of (j){Z, m, s) for 0 < m < K can be computed by the following equation 


K-l 

4>{Z, m, s) = 'f’iK, k, s) Tm{s) Um{Z, k — m,s) 

k—m 

+ ifiK, K + l,s) T(s) C/(Z, K-m,s) , (35) 

where the matrices UM{Z,k, s) and U{Z,k,s) are defined as 

UM{Z,k,s) = S{Z,{T^{s)S{Z,I,Tm{s)))\Tm{s)), 

U{Z,k,s) = SiZ,{TA{s)S{Z,I,TMis))f,T{s)). 

Proof. Using o and (1211), it follows that equation (1551) holds for m = if, where it is assumed that the value 
of the sum is zero. We prove by induction that it also holds for all m < K. We first derive a recursive equation 
satisfied by 4){-, m, s) in terms of (/)(■, m + 1, s). 

By substituting Tm{s) Z Tf^{s) for Z in (1551) we get an expression for 4>{Tm{s) Z Tf^{s), m, s), and subse¬ 
quently substituting this expression in the right-hand side of (I29L yields 

(j){Z, TO, s) ='ip{K, TO, s) Tm{s) + ipiK, TO, s) This) Z 
+ (s) Z Tj^^ (s), TO, s) (s) Z^ 

+ (j){TA{s) Tm{s) ZTffj^fs) T)2^(s), TO -I- 1, s) Ta{s) Tm{s) Z 
+ (fiTAis)ZTff\s),m + l,s)TA{s) (36) 

and iterating this equation leads to 


(j}{Z, TO, s) =ip{K, TO, s) Tm{s) ( Y ^m(s) Z^ I 

\h=0 ) 

CO 

+ ^ 0(Ta(s) Ttr{s)ZT^^{s) Tff\s),m + 1, s) Ta{s) T^s) Z'^ , 
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which can be rewritten as 


(j){Z, m, s) =ip{K, m, s) Tm{s) S{Z, J, Tm) (37) 

oo 

+ ^ </>(rA(s) TUs)ZT^\s) T^\s), m + l,s) Ta{s) T^{s) Z^ . 
h—0 

The recursive equation dSZl) is valid for m = 0,1,..., iiT — 1. 

We conjecture that for all m = 0,1,..., 77, the generating function fpiZ, rn, s) has the form 

K-l 

cl>{Z, m,s) = Y, fc, s) Tm{s) S{Z, y'=-™(s), Tm(s)) (38) 

k—m 

+ 77 + 1, s) T{s) SiZ, (s), T{s)) , 

so that (ESD follows by showing that the right expression for Y{s) is given by 

Yis)=TAis)SiZ,I,TM{s)) . (39) 

This conjecture will be proved by induction. We have already shown that it holds for m = 77. Now assume 
that it is valid for m + 1. To establish ((Mil for m, it suffices to prove, by virtue of dSil), that 

OO 

Y, HTa{s) TUs) Z T^\s) T^\s), m + 1, s)Ta{s) Z'^ 

/i^O 

K-l 

= Y i’{K,k,s)TM{s)S{Z,Y'^-^{s),TM{s)) (40) 

+ V;(77,77 + 1, s) T{s) S{Z, r^-™(s), T{s)) . 

It follows from Lemma [?731 that 

OO 

V-(77,77 + 1) T ^^(Ta Z T^\Y^-"^-\T) Ta Z^ 

h—0 

= ^(77,77 + 1) TS{Z, Ta S{Z, I, Tm), T) , (41) 

where we suppressed the dependence on s. Application of Lemma Id.81 is justified, since it is readily verified 
that the matrices in the above infinite sum satisfy the conditions mentioned in this lemma. Accordingly, for 
A: = m+1,...,77 — 1, and again suppressing the dependence on s, 

OO 

V'(77, k) Tm ^5(Ta tIj Z T^^ T]^\Y^-'^-\Tm) Ta T^ Z’^ 

h—Q 

= ^{K, k) Tm S{Z, yfc—1 Ta S{Z, 7, Tm), Tm) ■ (42) 

Combining (HTl) and we conclude, by virtue of the induction hypothesis, that (HUl) holds whenever y(s) 
satisfies (IMll . which completes the proof. □ 

Remark. Also in this case, as was already mentioned in Remark [21 the Laplace transform of the sojourn time is 
rational. This admits application of classical inversion techniques, yielding an explicit expression for the sojourn 
time distribution. In Section 13.21 we give an example of how to compute the density function of the sojourn 
time for a system with 77 = 2. 

3.1 Erlang inspection times 

In section [3] we assumed exponential inter-inspection times. In principle this can be extended to the case of 
phase-type distributed inter-inspection times [T], paying a cost in terms of model complexity. Indeed, in this 
case one should keep track, not only of the value of the service rate, but also of the phase of the inspection-clock. 
This translates into more complicated matrix expressions, but the basic logic of the computation of the sojourn 
time distribution remains the same. In fact, this is exactly the power of the proposed matrix generating function 
technique. For the sake of clarity and conciseness we are not going to treat here this extension in detail, but 
give a quick view of how it can be handled. 

We assume that the inspection times are Erlang(2,7) distributed. To keep trace of this we consider four 
states in the description of the system, {00,01,10,11}, where the first number specifies the speed of the system 
and the second the phase of the inspection clock. 
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Figure 4: Transition diagram for Erlang-2 inspection times 


The column vector 7r„ = (TToon, ■ ■ ■, satisfies (11811 with the following matrices 


M = 


Mo 0 

0 fii 


1 0 
0 1 


A = 


A 0 
0 A 


1 0 
0 1 


and Hi = A -|- r2, H 2 = M -\- K + T 2 and H 3 = M + A + T 3 , where 

r2 = 

The conditional sojourn times satisfy the following equation, see the corresponding formula 

a: 


/ 

7 

-7 

0 

-7 \ 


/ 

7 

0 

0 

0 \ 


-7 

0 

7 

0 

0 

7 

0 

0 

, Ta = 


-7 

0 

7 

-7 

0 

7 

0 

-7 

V 

0 

0 

-7 

7 J 


1 

0 

0 

-7 

7 / 




A -f Mi + 7 


Stj{n-l,m) w.p. Mi/(A-I-Mi + 7) 
5*^(71,771-1-1) w.p. A/(A-I-M* + 7) 
Sh{ij) {n, m) w.p. 7/(A -f Mi + 7) 


(43) 


with h{i,j) = h(i, j;n,m) = ((1 — j) ■ i + j ■ l{n + m > K},{1 — j)). It follows that the row vector 
(■000(’T-j TO; s),'ijjoi(ri, m, s), 'ipio(n, m, s), 0ii(n, m, s)) satisfies equation (1251) with the matrices H(s) = (s + 'j) 1 + 
A + M and 



/ 0 

7 

0 

7 \ 


/ 0 

0 

0 


ro = 

7 

0 

0 

0 

0 

0 

0 

0 

, ri = 

7 

0 

0 

7 

0 

0 

0 

7 


1 0 

0 

7 

0 ^ 


1 0 

0 

7 

0 J 


Since the matrix equations for the Erlang inspection times are similar to the exponential inspection times, 
all the subsequent matrix analysis in Proposition 13.21 and Theorem 13.41 are still valid. The value of the matrix 
i? is now given by 


i? = 


f Rii 

0 

0 

0 

R 21 

Rii 

0 

0 

R 31 

R 32 

R 33 

^34 

\ -^41 

Ra2 

R 3 A 

R 33 
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with 


i?ii 


R33 

R 34 


R 32 


Rai 


i?42 


7 + A + /4o - v^- 4 A/io + (7 + A + /ip)^ __ 7 -Rii _ . 

2 fj,Q ’ 7 + A + ^0 ~ ‘^fJ-oRii 

/xi (27 + 3 A + / 4 i) - (472 + (A - / 4 i )2 + 47 (A + ^1)) 

4^2 ; 

A _ 7 (i ?41 + ^21) + Ml (^ 32^21 + -^41^34) 

7 + A-m(-l + Kn+fl„) ■ 

7j?ii(-7 - A + Mi (-1 + i?ii + R33)) . 

-(7 + A - Mi (-1 + ^11 + R33W + (7 + ^1-^34)^ 

7-^21 (7 + Mi^ 34 ) (A^ — 2 Ami (—1 + R33) — Ml “ (~1 + -^33)^ + ^34)) 
(27 + A — Mi(~l + ^11 + R33 ~ ^34))^(A — Mi(~l + ^11 + R33 + -^34))^ 
_ 7-^21 (7 + M1-R34) (27(A - Mi (-1 + R33 + -^34))) _ . 

(27 + A — Mi(~l + ^11 + R33 ~ ^34))^(A — Mi(~l + ^11 + R33 + -^34))^ 

_ 7 -Rii (7 + Ml ^34)_ 

(27 + A — Mi(~l + -Rii + R33 ~ ^34))(A — Mi(~l + ^11 + R33 + ^34)) 


3.2 Analitical example 

In this section we briefly show that by using Theorem l3.41 we can get explicit expressions for the density function 
of the sojourn time in the system with inspection times, as highlighted in Remark [3l 

The computations are simple, but tedious as they require extensive use of matrix calculus, and usually it is 
easy to be assisted by symbolic computational software as we do for this example. 

To make computations easy, we wisely select the values of the parameters of the system such that all 
coefficients turn out to be rational. 

The parameters of the queue are 


Mo = 1; Ml = 3/2; 7 = 1/8; A = 9/8 . 


For the moment we do not hx the threshold, later we consider explicitly the case K = 2. The above choice of 
the parameters gives Rqq = 3/4 in (I^D) . The matrix R and the matrix function T{s) = M {H{s) — Fi — A)“^, 
are given by 

/ 3 0 A / —^ 0 

R={ t 3 )-, T{s)=( ^+ 3 ^- J_ 

V 4 4 / V (3+2s)(9+8s) 3+2s 

and the matrix functions T'a(s) = A {H{s) — Fi)”^ and Tm{s) = M {H{s) — Fi)”^ are equal to 


Ta{s) 


2 Tm{s) 

2(9+4s)(21+8s) 21+Ss J 


4 

9+4s 

_ 6 _ 

(9+4s)(21+8s) 


0 

12 


21+8s 


Solving the matrix system dSH) we get the following expression for S{R,I,Tm{s)), 


S{R,I,Tm{s)) 


9+4s n \ 

2(3+2s) ^ \ 

3 ( 3 + 5 ) 21+85 I 

2(3+25)2 4(3+25) / 


that allows the computation of the values of U{R,k,s) and UM{R,k,s) for any k > 0. As example we show 
such matrix functions for fc = 2, 


U{R,2,s) 


81(9+85) p, 

16(3+25)2(3+85) ^ 

81(69+88s) 81 

16(3+25)2(3+85)2 4(3+25)(3+85) 


Um{R, 2, s) = 


/ 81(9+45) 

( 32(3+2s)3 

81(30+115) 
\ 32 ( 3 + 25)^1 


» j 

81 ( 21 + 85 ) I 

64(3+25)2 / 


Remark. The expressions of S{R, I,Tm{s))), U{R,k,s) and UM(.R,k, s) do not depend on iF, so they can be 
used for any value of the threshold. The values of ■i/'(s), i/(n, 0, s), 7r„ and (/(i?, 0, s) in (1551) do depend on K via 
the respective formulas (1551) . (1551) . (155|) and (1551) . 

From here on we fix AT = 2. We have t^k = (3807/60644,1701/30322)3^ and after recursively computing 
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ip{k,0,s), for k = 1,2, we finally get V'(s), 


308367 13923657 44764461 130808703 

~ ~ 379025(3+ 2s)4 “ 9475625(3 + 2s)3 ~ 47378125(3 + 2s)2 “ 236890625(3 + 2s) 

14013 , 1587762 4755267 28797784929 

~ 15161(9 +4s) 9475625(11+ 4s)3 ~ 24636625(11 + 4s)2 “ 40034515625(11 + 4s) 

81216 ^ 18144 ^ 102060 ^ 55081053 

15161(3+ 8s)2 15161(3 + 8 s) 15161(9+ 8s)2 20497672(9 + 8 s) 

24064452 199526994 ^ 2950774277 , 90111 

“ 9475625(17+ 8s)3 “ 47378125(17 + 8s)2 1895125000(17+8s) 60644(21 + 8 s) 

the inverse-transform of which results into the following density function 

90111e-2“/8 14013e-®*/4 27e-3*/8(84 + 47f) 729e-9‘/®(75557 + 236601) 

■' ~ 485152 60644 15161 163981376 

+ 2436“^“/"^ (-1896150448- 1271985001+ 13803075i2)/2562209000000 

- (-11803097108 + 39905398801 + 1504028251^)/60644000000 

- 3e-®‘/2 (697646416 + 5968594801 + 2320609501^ + 214143751®) /7580500000 . 

Figure [S] plots the density functions of the sojourn time for K = 0,1, 2,3 using their exact expressions, instead 
of using numeric inverse transform as done later on in Section 21 



Figure 5: ^0 = 1; Mi = 3/2; 7 = 1 / 8 ; A = 9/8; and iF € {0, 1, 2, 3} 


3.3 First moment calculation 

Like in the previous section, we define u = E[S'] and i'n,m = E[S'(n, m)]. By taking derivatives in (l?51) and then 
computing the limit for s —> 0 we get 

^n,Tn (-^( 0 ) r l^n+m>K}) ~ — AL + 1^71,771+1 A + 6 , ( 44 ) 

where we used that H'{0) is the identity matrix. The vector e is the all-one vector. From (l26ll and after taking 
derivatives, we obtain 

n 

= e Y^{T{Q)f M-® (T(0))”-'=+® as n > 0 and m > iF , (45) 

fc=i 

with r(0) = M{H{Q) - Fi - A)"®, (r-i)'(O) = M"® and T'(0) = -T{Q)M-^T{Q). Here we used that the 
derivative of a matrix is given by 

n 

(H-”)' = A-^ . 

k^l 
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By Proposition 13.21 we can conclude that 


K-l 

i/ = ^ 0,0+) TT/f . (46) 

n—0 

From equation ([35]) we can compute 

K-l 

0, 0+) = ^ VK,k Tm(0) Um{Z, k, 0) + iyK,K+i T{0) U{Z, K, 0) 

k—0 

K-l 

-e Y, TM{0)U'M{Z,k,0)-eT{0)U'{Z,K,0) 

k^O 

K-l 

-e YTMiO)UM{Z,k,0)-eT'{0)U{Z,K,0) , (47) 

k=0 


with rM(0) = M{H{0) - Pi)-! and Tj^{0) = Tm( 0). The values U'j^{Z,k,0) and U'{Z,k,0) 

appearing in (IHHll can be computed by solving the following linear systems, see Lemma FA. II in the appendix, 

U'j^iZ, k, O) - Tm(0) U'^iZ, k,0)Z- T'j^,{0) Um{Z, k,0)Z = A'(Z, k, 0) 

[/'(z, k, 0 ) - r(o) u'{z, k,Q)z- r'(o) u{z, k,o)z = a'{z, k, o) 


with A{Z, k, s) = (TAis) S{Z, I, ^M(s)))^ 


4 Numerical experiments 


In this section we show some numerical examples, where we compute the stationary sojourn time distribution 
for a system with slow rate /zq = 1 and high rate fii = 3/2. 

In Figures [5] and 0 it is shown how the sojourn time distribution depends on the threshold K for the case 
of immediate switching times. In the first example, A < < Mii which implies that the system is stable for 

both service rates. Therefore, when K —>■ oo, one can appreciate that the sojourn time distribution approaches 
the one of an M/M/1 system with fixed service rate /zq (shown as dashed black line in FigureFB]). In the second 
example, we have A S [/zq, /zi). In particular, we have chosen A = /zq = 1, implying that the system approaches 
instability as AT —5> oo. 




Figure 7: A = I 


Figures IRl-im depict the sojourn time distribution for the case of exponential distributed inspection times. 
These figures refer to the case when A € /zi), and again, one can notice that as K ^ oo, the system becomes 

unstable. It is worth to notice that, when K = 0, the curve does not coincide with the M/M/1 with constant 
service rate /zi (shown as dashed blue line), since in the system with exponential switching, when inspection 
finds the system empty, the server switches to the slow rate and does not switch back till another inspection 
occurs. When 7 = 1000, the system switches almost immediately and therefore the sojourn time distribution is 
very close to the one of the pure M/M/1 system. 

In Figures [12] and [131 we plot again the results for A = 1, /zq = 1 and /zi = 3/2, but compare different 
values of 7’s. One can see that for 7 approaching A, the system behaves very close to a system with immediate 
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switching (shown as dashed black line). Indeed, for values of 7 > 1, one cannot distinguish the curve from the 
limiting one. This suggests that, checking the state of the system at a rate comparable to the arrival rate, can 
be considered from the point of view of the sojourn time as an immediate switching. This could be used in 
the design phase of the system, when balancing between costs (by reducing service rate) and performance (by 
increasing the service and inspection rate). 



Sojourn time distribution for: A = 1, /xq = 1, /ri 



Figure 13: K = Z 

= 3/2 and 7 e {1/100,1/10,1,10,100,1000, 00} 


5 Conclusions 

In this paper we studied the sojourn time distribution in an exponential single-server queueing system. Service is 
in order of arrival, and it is provided at low or high rate, which can be adapted at exponential inspection times, 
depending on the number of customers in the system. To determine the Laplace transform of the stationary 
sojourn time distribution, we proposed a new methodological tool, that is matrix generating functions. We used 
this tool to compute the Laplace transform of the sojourn time distribution in the system with inspection times. 
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Its expression is obtained recursively and shows a rational form that allows an immediate inverse-transformation. 
Numerical computations have shown, as expected, that if the inspection rate is large, the sojourn time of the 
system with inspections converges to the one of the system with immediate switching. 

We believe that the power of the matrix generating functions lies in its flexibility to analyze generalizations 
to phase-type services and inspection times. An interesting and promising direction for future research is to 
explore the applicability of this tool to analyze the more general class of quasi-birth-and-death processes [14] . 


References 

[1] S. Asmussen, Applied Probability and Queues, Springer, New York, (2003). 

[2] R. Bekker, S.C. Borst, O.J. Boxma, O. Kella, Queues with workload-dependent arrival and service rates, 
Queueing Systems, 46, 537-556 (2004). 

[3] R. Bekker, O.J. Boxma, An M/G/1 queue with adaptable service speed. Stochastic Models, 23, 373-396 
(2007). 

[4] R. Bekker, O.J. Boxma, J.A.C. Resing, Queues with service speed adaptations, Statistica Neerlandica, 62, 
441-457 (2008). 

[5] O.J. Boxma, B.H.B. Jonsson, J.A.C. Resing, V. Shneer, An alternating risk reserve process - Part II, Markov 
Processes and Related Fields, 16, 425-446 (2010). 

[6] J.W. Cohen, On the optimal switching level for an M/G/1 queueing system, Stochastic Processes and Their 
Applications, 4, 297-316 (1976). 

[7] J.W. Cohen, The Single Server Queue, North-Holland, Amsterdam (1982). 

[8] G.I. Falin, Aggregate arrival of customers in a one-line system with repeated calls, J. Ukrain. Math., 28, 
561-565 (1976). 

[9] G.I. Falin, On the Waiting-Time Process in a Single-Line Queue with Repeated Calls, J. Applied Probability, 
23(1), 185-192 (1986). 

[10] G.I. Falin, A Survey of Retrial Queues, Queueing Syst. Theory Appl., 7(2), 127-167 (1990). 

[11] D.P. Gaver, R.G. Miller, Limiting distributions for some storage problems. In; Studies in Applied Proba¬ 
bility and Management Science, 110-126 (1962). 

[12] J.M. Harrison, S.I. Resnick, The stationary distribution and first exit probabilities of a storage process 
with general release rule. Mathematics of Operations Research, 1, 347-358 (1976). 

[13] J. Keilson, L.D. Servi, A distributional form of Little’s Law, Operations Research Letters, 7, 223-227 (1988). 

[14] G. Latouche, V. Ramaswami, Introduction to Matrix Analytic Methods in Stochastic Modeling, SIAM 
(1999). 

[15] 1. Mitrani, Managing performance and power consumption in a server farm, Annals of Operations Research, 
202, 121-134 (2013). 

[16] M.F. Neuts, Matrix-geometric solutions in stochastic models - an algorithmic approach, Dover Publications, 
(1994). 

[17] R.W. Wolff, Poisson arrivals see time averages, Operations Research, 30, 223-231 (1982). 

A Technical proofs 

of Lemma \2.1\ . We can rewrite the expression ([5]) in the following form 

'0(n, m, s) = as{n + m) ip(n — 1, m, s) + hs(n J- m) ipin, m + l,s) . (48) 

With m = K — 1, equation © becomes 

^/(n, K — l,s) = Bsin + K — 1,1) tp{n, K, s) 

K-l 

+ as{n + k)Bs{n + K — 1, k — K + 1) ip{n — 1, k, s) 

k=K-l 

= bs{n + K — 1) K, s) + as{n K — 1) '0(n — 1, K — 1, s) 
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and therefore it holds true. Now assume that equation o is valid for m + 1. Then by (1481) , 

'0(n, m, s) =as{n + m) -ipln — 1, m, s) + bs{n + m) Bs(n + m + l,i4r—1 — m) K, s) 

K-l 


+ bs(n + m) as(n + k) Bs{n + m + l,k — m — 1 )'ip{n — l,k, s) 

k=m-\-l 

=as{n + m) Bs(n + m, 0) ip{n — 1, m, s) + Bs{n + m, K — m) tpin, K, s) 

K-l 

+ as{n + k) Bs{n + m,k — m)'ip{n — l,k,s) 


k=m-\-l 


K-l 

=Bs{n + m, K — m) K, s) + as{n + k) Bs{n + m,k — m) il){n — 1, k, s) . 

k—m 

where we have used the fact that the definition of Bs{k, h) implies that 

bs{k) Bs{k + 1, /i) = Bs{k, h + 1) . 


□ 


of Lemma \S.1[ By substituting in (l32ll the expression for S given in eiD we get 


BSZ = = '^B’^AZ^ 


h—0 


h—0 


-A=S-A 


which implies that the matrix 5" is a solution of the matrix equation. 

By assuming that S and S' are two solutions of this matrix equation, we would have that Y = S — S' \s the 
solution of the following system 

Y = ZY B. 


Iterating the last equation we get that 

Y = Z''YB'' n > 0 . 

This term converges to 0 as n —>■ cjo by the assumptions on the eigenvalues of the matrices Z and B. It follows 
that y = 0 and hence S is unique. □ 


of Lemma \,y.,yi The result follows from the following algebraic manipulations 

OO 

S{T 2 Z TfA, B) T 2 Z^ 

h—0 

OO OO 

= EE 5'= A (T 2 Z Tf ^ TfO'" T 2 Z^ 

h—0 k—0 

OO OO 

= EE 5'= A T 2 T’f Tf ^ Tff 1 T 2 T'f Z’^ 

h—0 k—0 

OO OO OO / OO \ 

= ^^2 Z^+^ = Yb'" AT2iYTi z''] 

^^0 k—0 k—0 \h—0 / 

OO 

= Yb^AT 2 S{Z, I, TOZ^ = S{Z, AT 2 SiZ, I, Ti), B) 

fc =0 

□ 


Lemma A.l. Let S{s) = S{Z, A{s), B(s)), then its derivative in s can be computed as the solution of the 
following linear system. 

S'{s) - Bis) S'{s) Z - B'is) Sis) Z = y4'(s) (49) 

Proof. By (15^ we have that 

SiZ, A(s + h), Bis + h)) — Bis + h) SiZ, A(s + h). Bis + h)) Z = A(s + h) (50) 

SiZ, Ais),Bis)) - Bis) SiZ, Ais),Bis)) Z = A(s) . (51) 

Subtracting the expressions above, adding and removing Bis + h) SiZ, A(s), Bis)) Z we have 
AS'(s) - Bis + h) AS'(s) Z - AB(s) 5'(Z, A(s), Bis)) Z = AA(s) 

with AS'(s) = SiZ,Ais + h),Bis + h)) — SiZ,Ais),Bis)) and similar notations for AA(s) and Ai?(s). Dividing 
for h and letting —>■ 0 the result follows. □ 
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